NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 


THESIS 


Approved  for  public  release;  distribution  is  unlimited. 


DTIC  QUALITS’  UrSPECTEDlT 


19981211  004 


REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction,  searching 
existing  <kta  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this 
burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services, 
Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management 
and  Budget,  Paperwork  Reduction  Project  (0704-0 1 88)  Washington  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

September  1998  Master’s  Thesis 

4.  TITLE  AND  SUBITTLE  BIAS  EFFECTS  ON  MOTION  STABILITY  OF. 
SUBMERSIBLE  VEHICLES 

5,  FUNDING  NUMBERS 

6.  AUTHOR(S)  Keith  L.  Payne 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 

Monterey  CA  93943-5000 

8.  PERFORMING 

ORGANIZATION 

REPORT  NUMBER 

9-  SPONSORING/MONTTORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSORING/MONTTORING 
AGENCY  REPORT  NUMBER 

1 1 .  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect 
the  ofScial  policy  or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 

12a.  DISTRBUnON/AVAILABILrrY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE* 

13.  (maximum  200  words) 

This  thesis  analyzes  the  nonlinear  characteristics  of  motion  stability  of  a  submersible  vehicle  in  combined  sway,  yaw,  and  roll  motions 
Previous  results,  at  zero  pitch  angles,  indicate  that  limit  cycles  are  generated  as  a  result  of  loss  of  stability.  In  this  work,  these  results 
are  extended  to  include  nonzero  pitch  angles.  This  analysis  can  determine  how  changes  in  vehicle  parameters  and  loading  conditions 
will  affect  its  operation  and  performance.  Stability  domains  are  generated  for  a  variety  of  vehicle  and  environmental  parameters.  A 
nonlinear  analysis  is  conducted  in  order  to  assess  the  stability  characteristics  of  the  resulting  limit  cycles.  The  results  can  lead  to 
design  guidelines  for  improving  vehicle  operational  envelopes. 

14v  SUBJECT  TERMS:  SurfaceAJnder  surface  vehicles.  Modeling  and  simulation. 

15..  NUMBER  OF 

PAGES  66  ■ 

‘ 

16.  PRICE  CODE 

17.  SECURITY  CLASSFICA- 
TION  OF  REPORT 

Unclassified 

18.  SECURITY  CLASSIFI¬ 
CATION  OF  THIS  PAGE 

Unclassified 

19.  SECURITY  CLASSIFICA- 
DON  OF  ABSTRACT 

Unclassified 

20.  LIMITATION  OF 
ABSTRACT 

UL 

NSN  7540-01-280-5500  .  Standard  Form  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std.  239-18  798-102 


1 


11 


Approved  for  pubKc  release;  distribution  is  unlimited 


BIAS  EFFECTS  ON  MOTION  STABILITY  OF  SUBMERSIBLE  VEHICLES 


Keith  L.  Pa)aie 

Lieutenant,  United  States  Navy 
B.S.,  Maine  Maritime  Academy,  1991 

Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 

MASTER  OF  SCIENCE  IN  MECHANICAL  ENGINEERING 

from  the 


NAVAL  POSTGRADUATE  SCHOOL 
September  1998 


ABSTRACT 


This  thesis  analyzes  the  nonlinear  characteristics  of  motion  stability  of  a 
submersible  vehicle  in  combined  sway,  yaw,  and  roll  motions.  Previous  results,  at  zero 
pitch  angles,  indicate  that  limit  cycles  are  generated  as  a  result  of  loss  of  stability.  In  this 
work,  these  results  are  extended  to  include  nonzero  pitch  angles.  This  analysis  can 
determine  how  changes  in  vehicle  parameters  and  loading  conditions  will  affect  its 
operation  and  performance.  Stability  domains  are  generated  for  a  variety  of  vehicle  and 
environmental  parameters.  A  nonlinear  analysis  is  conducted  in  order  to  assess  the 
stability  characteristics  of  the  resulting  limit  cycles.  The  results  can  lead  to  design 
guidelines  for  improving  vehicle  operational  envelopes. 
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I.  INTRODUCTION 


A.  DESCRIPTION  OF  THE  PROBLEM 

The  study  of  dynamic  responses  in  a  submersible  vehicle  using  a  nonlinear 
analysis  is  important  in  determining  operating  envelopes  for  the  vehicle.  Previous  studies 
have  shown  that  in  straight  line  motion,  the  coupling  of  the  sway,  yaw,  and  roll  equations 
produce  oscillating  losses  of  stability  in  the  system.  [Ref.  5]  Introducing  a  nonzero  pitch 
angle  to  the  vehicles  motion  will  allow  us  to  study  the  changes  to  the  stability  domain  for 
a  variety  of  environmental  conditions. 

In  our  work,  the  linearized  equations  of  motion  are  studied  using  an  eigenvalue 
analysis  to  determine  the  systems  stability  through  a  variety  of  operational  parameters.  A 
nonlinear  analysis  is  then  conducted  to  assess  the  stability  characteristics  of  the  resulting 
limit  cycles  and  their  impact  on  the  operating  characteristics  of  the  vehicle. 

B.  OBJECTIVES  AND  OUTLINE- 

In  this  thesis  we  expand  on  previous  thesis  work  which  examined  the  problem  of 
stability  in  straight  line  motions  of  a  submersible  vehicle  [Ref.  13].  The  primary  cause 
for  this  loss  of  stability  is  the  coupling  between  the  sway/yaw/roll  equations  of  motion  for 
a  submersible  vehicle.  We  know  that  the  loss  of  stability  creates  stable  limit  cycles  in 
straight  line  motion.  This  work  analyzes  the  effects  of  introducing  a  nonzero  pitch  angle 
to  the  generic  equations  of  motion  in.  order  to  determine  their  effect  on  the  creation  of 
limit  cycles. 
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The  model  used  for  this  work  is  a  variant  of  the  Swimmer  Delivery  Model  used  in 
[Ref.  5]  for  which  a  generic  set  of  hydrodynamic  and  geometric  properties  are  available. 


n.  EQUATIONS  OF  MOTION 


A.  COORDINATE  SYSTEM 

A  moving  coordinate  system  was  used  for  our  analysis  with  the  origin  fixed  at  the 
vehicles  center  of  buoyancy.  The  jc-axis  is  fixed  to  the  longitudinal  plane  of  symmetry  for 
the  vehicle,  the  y-axis  is  positive  to  starboard,  and  the  z-axis  is  positive  downward.  All 
symbols  used  in  the  development  of  the  equations  of  motion  are  summarized  in  Table  1. 

B.  GENERAL  FORM  OF  THE  EQUATIONS  OF  MOTION 

The  equations  of  motion  for  a  submersible  vehicle  in  the  horizontal  plane  are 
written  as  follows: 

Sway  equation: 

m[v  +  Ur-wp  +  Xc  {pq  +  r)-  (p^  +  )+  {qr  -  p)]= 

YpP  +  Yrr  +  Ypg  pa  +  Y^,qr  +  Y^Uv  +  Y^vw  + 

Y,U^6,+Y,v  +  Y^Up  +  Y,Ur  +  Y,^vq  +  Y,^wp  +  Y^,wr  + 

(W  -  fi)cos  0  sin  0  - 
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Yaw  equation: 


+  (/^  -  lyy  )pq  - (p"  -  )-  {pr  +  q)+  {qr  -  p)+ 

m[x(; {v  +  Ur-wp)-  -vr  +  w^)] = 

Npp  +  N^r  +  Np^  pq  +  N^^qr  +  NJUv  +  N^vm>  + 

NgU^S^  +  N^v  +  NpUp  +  N^Ur  +  N^^vq  +  N^wp  +  N^^wr  + 
[x^W  - XgB)cos0 sin 0  +  ();cW  -  yBB)sm 6  +  - 

P"'  [c^  h{xXv  +  xrf  +  C^^XxXw -  xqff^-^xdx 

K:,  >  '  UX^) 


Roll  equation: 

^xcB  +  (4  V  +  -4)- + 

mljyg  (w  - +  vp)- Zg  (v  +  f/r  -  w/?)]  = 

K^p  +  Kfr  +  Kp^pq  +  K^^qr  +  K^Uv  +  K^vw+ 

KpropU  ^+K,v  +  KpUp  +  K^Ur  +  K^^vq  +  K^wp  +  wr  + 

(>0^  -  >'gB)cos0  COS0  -  (zgW  -  ZgB)cosB  sin^ 

(3) 

The  rotational  velocity  equation  around  the  x-axis; 

(i>  =  P  (4) 

f/c/ denotes  the  cross-flow  velocity: 

Ucf  (^)  =  +  xrf  +  (w  -  xqf 
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Cn 

quadratic  drag  coefficient 

dr 

radder  deflection 

h(x) 

local  height  of  hull 

yy>^  zz) 

vehicle  mass  moments  of  inertia  about  body  axes 

xy  yz>^  zx  ) 

cross  products  of  inertia 

(K,M,N) 

moment  components  along  the  three  axes 

m 

vehicle  mass 

rotational  velocity  components  along  the  body  axes 

Euler  angles 

U 

constant  vehicle  speed  along  the  x  -axis 

{u,v,yv) 

translational  velocities  about  ix,y,z)  axes 

ix.y,z) 

distances  along  the  three  body  axes 

(XXZ) 

force  components  along  the  body  axes 

(XgJC’Zg) 

coordinates  of  the  center  of  gravity 

(XB,yB>ZB) 

coordinates  of  the  center  of  buoyancy 

fore  coordinate  of  vehicle  body 

^  tail 

aft  coordinate  of  vehicle  body 

,  Table  1:  Nomenclature 


C.  SIMPLIFICATIONS 

We  must  simplify  the  equations  of  motion  in  order  to  reflect  the  fact  that  we  are 
analyzing  motion  about  the  y-axis.  The  simplifications  that  we  employ  are: 


•  Acceleration,  w ,  in  the  z-direction  is  zero. 

•  ■  Acceleration  in  the  longitudinal  direction,  u ,  is  zero. 

•  Rotational  velocity,  q,  and  acceleration,  q ,  in  the  y-  direction  are  zero 
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D.  SIMPLIFIED  EQUATIONS  OF  MOTION 


After  applying  the  above  simplifications,  the  equations  of  motion  become: 

Sway  equation: 

ni^+ljr-wp+x^r-y^(p^  +r'‘)+Zoin= 

Y,p+Y,f+Yfiv^Y^vw+Y,V^d,+Y,v+y,Up*Y,Vr*Y^wp+Y„wr* 

(W -B)cos  Osin  (j)- 

h{xXv  +  xrf  +  Ci,  b{x)w^  (6) 

Ki,  ■  ^  ‘  U^[x) 


Yaw  equation: 

Izz^-IxyP"" -IyzPr  +  IxzP  + 
m[xc  (v + C/r  -  wp)- Jo  (f/ -  vr)]= 

MpP  +  N^r  +  N^Uv+N^vw+NgU^S^  +N.v  +  NpUp  +  N^Ur  +  N^wp  + 
wr + ,(xcW - XgB)cos 0 sin 0  +  J b 5)sin $  +  - 

1'"°”  [c^  ;i(j:Xv  +  xrf  +  Cp.  b{x)w^ 

>  ‘  U^[x) 


Roll  equation: 

K^p  +  K,f+Kfiv  +  K^vw+  K^JJ^+K,v+  K^Up  +  K,Ur  + 
^wp'^P  +  +  (3^0^  ~  ^  ^  ~  ZbB)cos  6  sin  0 
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Roll  rate: 
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ni.  LINEAR  ANALYSIS 

A.  LINEARIZATION 

The  simplified  equations  of  motion  can  be  written  in  the  matrix  form: 

Ax  =  Bx+g(x)  (10) 

where  the  state  vector,  x  ,  is  defined  as, 

V 

r 

x  = 

P 

and  the  state  matrices  are, 


m-Y.  mXa-Y^  -mZc-Y^  0 

fnXc-N^  -I„-Np  0 

-mZa-K,  I„-K,  0 

0  0  0  1 

and 


Y,+Y^w  Y^U+Y^,-mU  Y^U  +  Y^w+mw  O' 

N,+N^w  -mXaU  +  N^U  +  N„^w  N^U  +  N^w  0 

KJU  +  K^w  mZcU +  K^U  ■>rK^^w  -mzGW+ K^U  +  K^w  0 

0  0  10 
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The  g(x)  term  remains  as  given  [Ref.  13]. 

The  nonlinear  terms  are  linearized  around  a  nominal  point: 

Xo  =[vo,^o>Po’0or  =0 


A  Taylor  series  expansion  is  applied  to  the  nonlinear  terms  about  the  nominal  point,  Xo , 
and  keeping  only  the  linear  components,  the  equations  of  motion,  written  in  matrix  form, 
become: 


A  x  =  B  X 


(11) 


where. 


and 


A'  =  A 

YU 

Yp-mU 

N,U 

—mx^U  +  N^U 

K^U 

mZaU  +  Kp 

0 

0 

YpU  (W-B)cos& 

NpU  {x^W -XsB)cose 
KpU  {-ZGW  +  ZBB)cose 
1  0 


To  assess  the  dynamic  stability  of  the  vehicle,  an  Eigenvalue  analysis  is  performed  in  the 
next  section. 
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B.  LOSS  OF  STABILITY 


An  Eigenvalue  analysis  is  used  to  determine  the  stability  of  the  linearized  system. 
Stability  is  dependent  on  the  location  of  the  Eigenvalues  or  the  roots  of  the  characteristic 
equation: 

det(jB'-M'  =  0)  (12) 

in  the  polynomial  form: 

AA^+BA^+CA^+Z>A  +  £:  =  0  (13) 

The  coefficients  equation  (13)  are  given  in  [Ref.  13,  pg.  1 1].  Using  Routh’s  criterion  we 
can  examine  the  stability  of  the  system.  The  following  conditions  must  be  applied  to  the 
characteristic  equation  (13)  to  ensure  all  roots  have  negative  real  parts: 

BCD-AD^-EB^  >0  (14) 

E>0  (15) 

If  E  is  less  than  zero,  one  real  root  of  (13)  becomes  positive  and  the  system  will  become 
unstable  in  a  divergent  manner  [Ref.  9].  This  is  the  case  of  a  directionally  unstable  ship 
which  is  well  known  in  the  literature  [Ref.  3].  If  the  condition  in  (14)  is  not  met,  it  means 
that  there  is  at  least  one  complex  conjugate  root  with  real  parts  and  will  result  in  an 
oscillatory  response,  indicating  loss  of  stability. 
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To  analyze  the  limiting  case  of  loss  of  stability  equation  (14)  must  be  solved  in  the 


following  form: 


BCD-AD^  -EB^  =0 


(16) 


The  result  of  this  equation  will  be  a  curve  of  zg  as  a  function  of  xg  and  will  be  our  locus 
for  loss  of  stability.  We  can  express  the  coefficients  of  equation  (16)  in  the  form: 


A  —  A^Zq  “1“  A2ZQ  *1”  A^ 

(17) 

B  =  B^Zq  +  B2ZQ  +  B^ 

(18) 

A\,  A2,  and  A3  are  of  the  form  given  in  [Ref.  13,  pg.  14]. 

B\  and  B3  are  of  the  form  given  in  [Ref.  13,  pg.  14].  .With  the  addition  of  a  pitch  angle,  w. 


B2  takes  the  form: 

S,  =-m(Ar,J/X/„ -W,Xf,C7) 

+  mYp{UN^  -UmXc)+mK,{UN^  -Unvc^ ) 

-m(Ar.  +I„XUY,  -Um)+m{N -Y,) 

-  m{K,  - 1 ^  iN^U  )+  mC/F,  [mx^  -N,) 

- mU(N^  +I^\m- Y, )+  mUK, {mx^  -N,) 

+  mZow(m-F^X^a -^f  )-fnZcw{^G -^rX^c  “•^v) 

C7  =  CjZg  +  C2Z(j  +  C3  (19) 
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Cl  and  C3  are  of  the  form  given  in  [Ref.  13,  pg.  15].  With  the  addition  of  a  pitch  angle, 
w,  C2  takes  the  form: 

C2  =mU{mXa -N,\Y^U)-mUYp{N,U)-mUKXN,U) 
+mU{N^+ljY,U)-mU{N^uXm-Y,) 
+W{m-Y,ll,-N,)-W{wxa-Yjmx^-N,) 

-m{XsB-XaWXmXa  -7.  )+m{UN^ -Unvc^  l^pU) 

-m{UY^  -UmiN pU)+m{KJj'pN ^  -Umx^) 

-mZaM>{rn-Y,\UN^  -Umx^  )-mZow{Y^U\l ^ -N,) 

-mzc'w{fnx(; -YfX^^U  )+mZcw(7nXc -N^  -Um) 

D  =  D,Zg+D^  (20) 

D\  is  of  the  form  given  in  [Ref.  13,  pg.  16].  With  the  addition  of  a  pitch  angle,  w,  D2 
takes  the  form: 

D,=UKX^,B-x„wXm-Y,)-UKXN.UtYpV] 

+UKM,vlr.U)-(K,-I,Xx,B-x„wiY,v) 

+  K,(xiB-x^W'iUY,-Vm)-(K,U'ix,B-XoW\mx^-Y,) 

+  (K,V\UN,  -UmXo1ffl)-(K,vyfrr, -VmtN ,u) 

-{k,uIy,uXUN, -Umx^)+(K^U\UK,  -UmlN,U) 

+mzgw(r^l7 XP^r  -UmX(^  )-mZgw{fJK^-UmXN JJ ) 

The  equation  for  the  coefficient,  E,  remains  unchanged  [Ref.  13,  pg.  16].  Applying  the 
stability  criterion,  equation  (16),  and  utilizing  them  in  the  resulting  fifth  order 
polynomial,  F,  [Ref.  13,  pg.  18]  we  are  able  to  solve  F  using  the  MATLAB  program  in 
Appendix  A.  The  curves  show  zq  as  a  function  of  xg  and  we  show  results  for  varying 
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pitch  angles,  w,  and  varying  forward  velocities,  U.  On  all  of  the  following  graphs  the 
pitch  angle  is  varied  from  10  degrees  to  -10  degrees  in  increments  of  five  degrees.  The 
top  curve  is  the  10  degree  curve  and  the  bottom  curve  is  the  -10  degree  curve. 


ZG  in  ft  ZG  in  ft 


Figure  3:  ZG  vs.  XG  for  U  =  5  ft/sec 


Figure  4:  ZG  vs.  XG  for  U  =  6.5  ft/sec 
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Figure  5:  ZG  vs.  XG  for  U  =  8  ft/sec 
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From  the  results  of  Figures  1  through  5,  the  following  conclusions  can  be  drawn: 

1.  In  all  cases,  a  sufficiently  high  metacentric  height  is  required  in  order  to  ensure 
vehicle  stability.  Regions  of  parameters  that  fall  below  the  critical  curves 
correspond  to  dynamic  instability. 

2.  The  critical  metacentric  height  that  is  required  for  dynamic  stability  is  an  increasing 
function  of  both  vehicle  speed  and  trim  angle. 

3.  Static  stability  alone  does  not  necessarily  ensure  dynamic  stability  of  motion  during 
the  turn. 

4.  The  loss  of  stability  experienced  here  is  a  dynamic  loss  of  stability.  At  the  critical 
metacentric  height,  one  pair  of  complex  conjugate  eigenvalues  possesses  a  zero  real 
part.  This  is  an  oscillatory  loss  of  stability,  which  can  not  be  predicted  by  considering 
the  uncoupled  sway/yaw  equations  of  motion  alone. 
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IV.  NONLINEAR  ANALYSIS 


A.  INTRODUCTION 

In  the  previous  chapter  we  performed  a  linear  analysis  of  the  equations  of  motion, 
observing  that  with  changes  to  specific  parameters  (a:g,  zq,  w)  it  is  possible  to  pass 
through  a  stable  region  to  a  region  of  loss  of  stability. 

In  previous  work  [Ref  13]  it  was  shown  that  these  bifurcations  to  periodic 
solutions  were  all  supercritical.  This  means  that  limit  cycles  were  produced  after  the  loss 
of  stability.  By  introducing  a  pitch  angle,  w,  in  the  nonlinear  analysis  we  will  analyze 
whether  the  bifurcations  to  periodic  solutions  will  remain  supercritical  and  what  changes 
occur  to  the  limit  cycles  themselves. 

B.  TfflRD  ORDER  EXPANSIONS 

Our  linearized  system  was  written  in  the  form  of  equation  (11),  ignoring  the 
nonlinear  terms.  Including  the  nonlinear  terms  changes  the  form  of  equation  (1 1)  to: 

A'x  =  B'x+g(x) 


where: 
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8  Ax) 
sAxX 


Keeping  terms  up  to  third  order,  the  matrix  g(x)  can  be  written  in  the  vector  form: 


^  ^  {x)+  8  {x)+  c(x) 


(22) 


where  g  (x)  contains  the  second  order  nonlinear  terms  and  g  (x)  contains  the  third 
order  nonlinear  terms.  The  cross  flow  integrals  and  the  second  order  nonlinear  terms 
remain  unchanged  with  the  addition  of  a  pitch  angle,  w.  [Ref.  13,  pg.22,  23]  However 
the  third  order  nonlinear  terms  take  the  form: 

where: 

gp)=-/f)-l(W-S)/»'cos0 
/  6 

gf^=-I»)-l{xW-x,B)li\ose 

o 

g?=^{z,W-z,Bycose 

«?'=0 

The  Taylor  series  expansion  yielding  the  second  and  third  order  linear  terms  of  the  cross 
flow  integrals,  and  the  inverse  of  the  system  matrix,  (A')"'  ,  remain  unchanged.  [Ref.  13, 


g(x)= 
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pg.  24, 25].  However  the  B’  matrix  has  changed  as  shown  in  the  Linear  Analysis  section 


of  this  work.  These  changes  to  the  4  x  4  F  matrix  are: 


Eil 

Ill. 

IlL 

IlL 

D 

D 

D 

D 

Fn 

F 

^  22 

F 

^23 

I2L 

D 

D 

D 

D 

F,: 

F32 

F33 

^34 

D 

D 

D 

D 

0 

0 

1 

0 

where: 

F„  =an{Y,U)+<.„{N.Uha,,{K,U) 

Pi2  =  “ii fr.  - '"X'  +<’1! - "“o + 

F„^a„{Y,u)+a,,(N,,u)+a„{K^u) 

Fj4  =  — B)cos9  +  ai2{x(;W — XgB)cos6  +  ai^{—  ZgW  +  ZBB)cosd 

F„=a,,{Y.U)+ajN,U)+a^iK,U) 

Fzi  =0!i(l',-'"X'+<'2i(W,-mxr.3)y+a„(X',+mz.;)[/ 

“‘n{Y,v)+ajN^v)+a,,(K„u) 

F24  =  ^21  (^  —  .S)cOS 0  +  022  0  +  023  (~  ZqW  +  Zg B)cOS  0 

Fj,  =‘‘AY.U)+‘‘„{N,V)+a„{K,U) 
Fi,=a,,{Y,-mp+a^,{N,-mx^'fJ+ay^(K,+mz^yj 
F„  =a„(Y,u)+ajN^u)+a„{Kfi) 

F34  =  O31  (W  -  B)cos  0  +  O32  W  -  X5B)cos  0  +  O33  (-  ZgW  +  Z5  B)cos  0 
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The  remaining  elements  of  the  nonlinear  and  constant  terms  remain  unchanged  [Ref .  13, 
pg.  28,  29]. 

C.  COORDINATE  TRANSFORMATIONS 

To  continue  our  analysis  it  is  necessary  to  bring  our  transform  our  coordinate 
system  from  state  space  to  a  normal  coordinate  system.  This  transformation  is  performed 
in  the  manner  given  [Ref.  13,  pg.  29,  30]. 

D.  CENTER  MANIFOLD  EXPANSIONS 

The  center  manifold  expressions  are  of  the  form  given  [Ref.  13,  pg.  30-33]. 
However,  the  coefficients  /y,i=l,  2, 3;  j=5, 6, 7  are: 

K5=^ya(^3i+^li) 

+/y,m3,m2,  +yc'”2i"^u) 

+  ygWml-ysBml^ 

ke  =^yo  (2"i3l"»32  +  2^2,  +  ^22  ) 

* 

+/y,(m3,/n22  +m32m2, ) 

+  yc(m2,m,2+»i22m,i) 

/^.(m3im22  +W32m2i 

+  myo  (m„m32  +mi2m3, )+  2{y(;W  -  5)^41  "*42 


D 
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ki  =  ^yG(^k+fnli) 

"*  '^{^xy^32  ^  yz^31^'n  ■*"  J  G^22^\l  ) 

-^[jxy^32^22  +  ^  yz^h  +'WyG«*12'«32  +  (^G^  "  ] 

T  ^21  /  2  2  \ 

4.5  =  — 3'GKl+'”2lj 

+  +^>^"^3l'«21  +>’G'”2l'«n) 

«23  f +^)-z'”2^  +W*yG'«n'”31  ^ 

^  [+}'G^"3«-yfi5'«4^  ; 

4.6  =  ^  >'g  (2«^3i'”32  •+  2wi2,  +  m22  ) 

fl„  P^;c>'”3l"»32+^yz("*3l"^22+'«32'«2l)' 

+  -^ 

L+>’c("^2i'«12+"*22"Iii) 
a23  r ('”3i"^22  +  "J32"«21  )+  yz"^2lfn22 
^  L+^3’G(^ll»^32+"*12"*2l)+2(}'G^->/r^Kl'”42. 

I  ^21  ^2  2  \ 

4.7  ~  3^G  1^32  ■*■  ^22  j 

+  ^(^^"^34  +  Iyz^3im22  +  y 0^21^12  ) 

-■ ^t;ty'”32'«22  +  ^ yz^22  +  "*>'g'”i2'”32  +  (^'g^  “  )’fi^)”42  ] ' 

4,5  ~  ^’g  (”^3l  ^21  ) 

+  ^  +  ^yz'W3:'W21  +  3'G'”2l'Wn  ) 


(25) 


(26) 


(27) 


(28) 
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(29) 


+/y,m2,  +mycm^^m3^  +ycWml^-yBBml^} 


D 


h,6  =-^3'g(2"^3i"^32  +2"l2I  +"^22) 


+  - 


a 
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(2 
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D 


2/^m3,m32  +/3.3(m3,m22  ) 

+  yc("J2i"ii2+'«22"»n) 

^;ty('”3l'”22  +"^32"J21  )+2/3^m2,m22 

+  me  i^nmi  +  '”l2'”3I  )+  ^{ya^  -  >'b^)^4i'”42  J 


1  ^31  /  2  2  ^ 

^3,7  ~  3^G  1^32  ^22  / 

+  ^(^^>"^3^  +'^yz"^32"^22 +>'g"J22"Ii2) 

-■ ^t;^"*32'«22  +^jj«i22  +"iyc"*12"*32  +  “  >’b^)^42  ] 


(30) 


(31) 


E.  AVERAGING 

The  procedure  for  averaging  the  equations  up  to  the  third  order  is  the  same  one 
given  in  [Ref.  13,  pg.  37-40].  The  addition  of  a  nonzero  pitch  angle  yields  new 
coefficients  /y,  7=1, 2, 3;  j=l,  2,  3, 4.  They  are; 
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7  —  ^11 

"  D 


C 


(jEqWZj^,  +3£jmf,m2i  +3£2^n"^2i  +-^3^21) 

+-(l¥-5)cos&n4^ 


a 


12 


D 


67 


(Ejmfi  +3£^2^n'”2i  +3£3m„m2i  +£47722,) 


+  —  (xgW  -  JCg  jB)cOS  ^41 


OI/ 


(32) 


a 


11 


M2 


D 


67 


3£o772,^j772,2  +3£,  (772,^22  +  2772,, 772,27722,  ) 

+  3£2  (772,27722,  +27722,77222772,,  )+ 3£3 7722,772- 


+ — (W  -  jB)cos  037724, 

m^2 


a 


12 


D 


Cj^  r3£,  772,,  772,2  2^p^\l^22  +  2772,, 772,27722,  )+ 3£j  (772,27722,  +27722,77222772,,  )')' 


67 


+  3E^ml\m22 


+  —  {xQW-XgB)  cos  037724,  ^42 


it-  (zg^  ~  Zft£)cos  d3ml^m 

6D 


'42 


(33) 


/  —  _.£il 
D 


r 3£o772„772,2  +3£,  (772,27722,  +  2m^^m^2^22  )+  '^^2  ^22'^n  +  '^^21^22^^12 


67 


+  3£3  7722,7724 


+  — (W -0)cOS  037724,772, 
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a 


12 


D 


— +3£2^n^2i  +  +  £4/7222  )h — {xgW  —  XgB^cosBm, 

67  6 


+7^fe<;M'-z,B)cosftn= 

bu 


42 


(34) 


7  -  ^21 


a 


22 


D 


1  -_^ 

D 


— ^(£o/n.^j  +3£,/nf,/n2,  +3£2m„m2,  H-Ej/Wj, )+-7(W-£)cos^4] 
67  6 


(Ej/w-Tij +3£2/72]j//i2i  ■*"3£3//Ii]^2i '^■^4^21  ^  X gS'^cosOiti^^  (35) 


67 


67 


3Ef,mfim^2  +3£,(mf,/n22  +2/?2i,/n,2/M2i)+3£2(/Wi2/W2i  +2/W2i/n22mj] ) 


^  + 3^3  7712177722 


+— (W -B)cos03/«4,/n, 
6 


42 


a 


22 


£> 


67 


3E^m^^m^2  +3£2(/wf,/n22  +2/ni]/ni2/M2i  )+3£3(/Wi2W2i  +  2/^21  "^22  "^11) 


+  3E^ml^m22 


+—{xqW -XgB)cos,d3ml^m^2 


+-^(zjW-ZB£)cos03/n4,/n, 


■42 


(36) 


1  —  —  ?2L 

23  -  ^ 


^23 

Uy 


67 


3£o/n,,/n,2  H-3£,  (/n,^2^2i  +2/Wn/Mi2Wi22)+3£2(/n22///ii  +2/n2,/n22/ni2) 


+  3£3/n2,/n^ 


+-(W-S)cose3/n4i  77^42 

6 


a 


22 

Z) 


Cq  f3£,/niimf2  +3£2(/nf2/n2,  +2/nii/nj2/n22)+3£3(/n22^ii  +2/n2,/n22///i2) 
^+3£4/n2,/n22 

+—{xqW — Xg£)cos03/M4,/n42 

6 


26 


27 


^33  ” 


a 


31 


D 


Cd,  f  +  3£:,  ^ 

+  3^2  (^22^11  ++2m2im22W,2  )+3£:3  m^im 


67 


2 

22  ) 


+  i(W-B)cos03m4,  m\ 

6 


a 
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D 


■Dr 


67 


1 


3£’,mi,mf2  +3£:2(m,2'”2i  +  2w„mi2m22) 
+  2m2im22m,2 )+ 


H — -  XgB)cos  dSm^^m 
6. 


+  -^tGl^-^B^)cos03m4im42 

oL' 


(41) 


/  —  ^31 

D 


-^{Eotnl^  +3£',mf,m2,  +3£:2W,2m2\  +£3m22) 


+  ^(W-B)cos0m 


“32 

D 


— —{Eyml^  +3£'2^h"^2i  +3^3mj2m22  +jE'4m22) 

67 

+  — (aTjjW'  -  j[:gB)cos  0m 42 
6 


oJD 


(42) 


The  remaining  procedure  for  determining  the  equation  for  the  radius  of  the  resulting  limit 
cycles  is  identical  to  the  one  described  in  [Ref.  13,  pg.  43, 44],  which  is: 


R  =  a'£R  +  KR^ 


(43) 
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F.  LIMIT  CYCLE  ANALYSIS 

At  steady  state,  ^  =  0 ,  equation  (43)  becomes: 

0  =  R{a'e  +  KR^) 

Equation  (66)  has  two  solutions: 


R  =  0 


(45) 

(46) 


Equation  (45)  is  the  trivial  solution  and  gives  no  useful  information.  Equation  (46)  gives 
a  constant  amplitude  limit  cycle,  R,  in  the  cartesian  coordinate  system.  This  limit  cycle 
will  exist  if  the  quotient  inside  the  radical  sign  is  positive: 

(47) 


K 


This  condition  is  necessary  for  to  be  a  real  number.  Since  a’  is  always  negative,  the 
existence  of  limit  cycles  depends  on  the  parameter  K.  The  introduction  of  a  nonzero  pitch 
angle  does  not  change  the  dependence  the  limit  cycle  has  on  the  parameter  K.  Stated,  this, 
dependence  is: 

•  If  K<0,  periodic  solutions  exist  and  they  are  stable. 

•  If  K>0,  periodic  solutions  exist  and  they  are  unstable. 
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G.  RESULTS  AND  DISCUSSION 


Results  for  the  stability  parameter,  K,  are  presented  in  Figures  6  through  15 
Figures  6  through  10  provide  results  for  K  at  a  constant  forward  speed,  V,  and  varying 
pitch  angles.  The  pitch  angles  were  varied  from  positive  10  degrees  to  negative  10 
degrees  in  5  degree  increments.  In  Figures  6  through  10,  the  bottom  curve  represents 
solutions  for  positive  10  degrees  and  the  top  curve  represents  negative  10  degrees.  It  is 
clear  that  all  values  of  K  are  negative  indicating  they  are  stable  solutions.  Notice  that 
decreasing  pitch  angles  the  solutions  for  ^become  less  negative,  indicating  that  while 
they  are  stable,  these  solutions  are  less  stable  than  those  at  the  higher  pitch  angles.  In 
Figures  1 1  through  15,  solutions  for  the  stability  parameter,  K,  are  again  represented,  but 
with  the  pitch  angle  held  constant  and  the  forward  speed,  U,  varied  from  2  ft/sec  to  8 
ft/sec  in  1.5  ft/sec  increments.  The  bottom  curve  in  Figures  1 1  through  15  represents  U  = 
2ft/sec.  As  forward  speed  increases,  the  curves  representing  the  stability  parameter  K 
tend  to  become  more  negative,'but  in  a  more  pronounced  fashion  than  those  where  U  was 
held  constant. 
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Figure  6;  XG  vs.  K*GAMMA  for  U=2  ft/sec 
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K'GAMMA  K'GAMMA 


U=5  ft/sec 


Figure  8:  XG  vs.  K*GAMMA  for  U=5  ft/sec 


U=6.5  ft/sec 


Figure  9:  XG  vs.  K*GAMMA  for  U=6.5  ft/sec 
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Figure  10:  XG  vs.  K*GAMMA  for  U=8  ft/sec 


33 


K'GAMMA  K'GAMMA 


THETA=10deg 


XG 

Figure  12:  XG  vs.  K*GAMMA  for  THETA=5  deg 
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K*GAMMA  K*GAMMA 


1 


THErA=0  deg 


xIO"® 


Figure  13;  XG  vs.  K*GAMMA  for  THETA=0  deg 


Figure  14:  XG  vs.  K*GAMMA  for  THETA=-5  deg 
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THETA=-10deg 


Figure  15:  XG  vs.  K*GAMMA  for  THETA=-10  deg 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  presented  a  continuing  study  of  the  formation  of  limit  cycles  due  to  the 
coupling  of  the  sway/yaw/roll  equations  of  motion.  We  have  shown  that  loss  of  stability 
occurs  in  the  form  of  stable  limit  cycles  and  that  the  addition  of  a  nonzero  pitch  angle  does 
not  significantly  affect  the  formation  of  these  limit  cycles.  Through  a  linear  analysis  of  the 
sway/yaw/roll  equations  of  motion  we  demonstrated  that  the  addition  of  a  nonzero  pitch 
angle  affects  the  domain  of  stability  of  straight  line  motion,  especially  at  higher  speeds. 
This  was  validated  by  the  nonlinear  analysis  as  well.  As  a  recommendation  for  further 
study  in  this  area  we  suggest  that  the  analysis  be  expanded  to  include  coupling  into  the 
heave/pitch  directions  of  motion  as  well  as  the  effects  automatic  path  control. 
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APPENDIX 


The  following  is  a  list  of  the  MATLAB  and  FORTRAN  programs  used  in  this  thesis. 
Complete  printouts  of  the  programs  accompany  this  list. 

•  THESIS1.M 

A  MATLAB  program  for  performing  the  linear  analysis  section  of  this  thesis. 

•  HOPF_lNEW.FOR 

A  FORTRAN  program  for  performing  the  nonlinear  analysis  section  of  this  thesis. 
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%  THESIS  1.M 

% 

%  LOSS  OF  STABILITY 

^  *****H!******=I=**************************** 

clear 

a=l; 

W=12000; 

IXX=1760;  IYY=9450; 

IZZ=10700;  rXZ=0;  IXY=0; 

IYZ=0;  D=17.425;  RHO=1.94; 

G=32.2;  U=8.0;  M=W/G;  B=W; 

THETA=0; 

THETA=THETA*pi/l  80; 

OMEGA=U*tan(THETA); 

ND1=0.5*RHO*L^2; 

%  DEFINE  HYDRODYNAMIC  COEFFICIENTS 

YPDOT=  L270e-04*ND  1  *L''2; 
YVDOT=-5.550e-02*NDl  *L; 

YRDOT=  1 .240e-03*ND  1  *L'^2; 
YP=3.055e-03*NDl*L; 

YPOMEGA=0; 

YP=YP*U+YPOMEGA*OMEGA+M*OMEGA; 
YV-  9.310e-02*NDl; 

YVOMEGA=0; 

YV=YV+YV0MEGA*0MEGA; 
YR=-5.940e-02*ND  1  *L; 

YROMEGA=0; 

YR=YR+YROMEGA*OMEGA; 

NPDOT=-3.370e-05*NDl  *L^3; 

NVDOT= 1 .240e-03  *ND  1  *L'^2; 
NRDOT=-3.400e-03*NDl*L''3; 

NP=-8 .405e-04*ND  1  *L''2; 

NPOMEGA=0; 

NP=NP+NPOMEGA*OMEGA; 

NV=- 1 .484e-02*ND  1  *L; 

NVOMEGA=0; 

NV=NV+NVOMEOA*OMEGA; 

NR=- 1 .640e-02*ND  1  *L^2; 

NROMEGA=0; 

NR=NR+NROMEGA*OMEGA; 
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KPD0T=-1 .01e-03*NDl  *L^3; 

KVD0T=1 .27e-04*ND  1  *L^2; 

KRDOT=-3.37e-05*NDl  *L'^3; 

KP=- 1 . 1 0e-02*ND  1  *L^2; 

KPOMEGA=0; 

KP=KP+KPOMEGA*OMEGA; 

KV=3.055e-03*NDl*L; 

KVOMEGA=0; 

KV=KV+KVOMEGA*OMEGA; 

KR=-8.41e-04*NDl*L''2; 

KROMEGA=0; 

KR=KR+KROMEGA*OMEGA; 

flag=0; 

forXG=0:0.01:2, 

flag=flag+l; 

xg(flag)=XG; 

a=IXX-KPDOT;  b=KP*U;  e=KV*U; 
f=KRDOT;  i=YP*U;  j=M-YVDOT;  k=YV=^U; 
l=XG*M-YRDOT;  m=U*(YR-M);  o=NPDOT; 
p=NP*U;  q=-XG*W;  i^XG*M-NVDOT; 
w=U*(NR-XG*M);  x=NV*U;  u=IZZ-NRDOT; 

al=-u*M^2; 

a2=-u*M*YPDOT-u*M*KVDOT+M*o*l+f*r*M; 

a3=a*j*u-a*l*r-u*KVDOT*YPDOT+KVDOT*o*I+f*r*YPDOT-f*o*j; 

bl=(w*(M'^2))+(r*U*(M'^2)); 

b2=-M*e*u-M*u*i+w*M*YPDOT+w*KVDOT*M'-o*m*M+p*l*M- 
f*x*M+r*M*U*YPDOT-o*j*M*U+r*KR*U*M; 
b3=-e*u*YPDOT+e*o*l-u*i*KVDOT+w*KVDOT*YPDOT- 
KVDOT*o*m+KVDOT*p*l-a*j  *  w-a*k*u+a*l*x+a*r*m-b*j  *u+b*l*r+f*r*i- 
f*x*YPDOT-K)*k*f-f*p*j+r*KR*U*YPDOT; 

b2=b2+M*OMEGA*j*u-M*OMEGA*l*r; 

cl=-x*(M^2)*U; 

c2=r*i*M*U-x*M*U*YPD0T-x*KR*U*M-K)*k*M*U-p*j*M*U+j*u*W-l*r*W- 
q*l*M+w*i*M-m*p*M+e*w*M;  • 

c3=a*k*w-a*m*x+b*j*w+b*k*u-b*l*x-b*r*m+r*i*KR*U- 

x*KR*U*YPDOT+o*k*KR*U-p*j*KR*U-q*l*KVDOT+w*i*KVDOT-m*p*KVDOT- 

e*u*i+e*w*YPDOT-e*o*m+e*p*l+f*q*j-f*x*i+f*p*k; 
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c2=c2-M*OMEGA*j*W-M*OMEGA*k*u-M*OMEGA*l*x+M*OMEGA*r*m; 


d2=q*j*M*U-x*i*M*U+p*k*M*U4q*m*M-j*w*W-k*u*W+l*x*W+r*m*W; 

d3=q*j*KR*U-x*i*KR*U+p*k*KR*U-f*q*k+q*m*KVD0T-e*q*l+e*w*i-e*m*p- 

b*k*w+b*m*x; 

d2=d2+M*OMEGA*k*w-M*OMEGA*U*(KR-M)*x; 


e2=k*w*W-m*x*W-q*k*M*U; 

e3=e*q*m-q*k*KR*U; 

f5=(-e2*(bl''2))+b  1  *cl  *d2; 

f4=((b2*cl+bl*c2)*d2)+bl*cl*d3-al*(d2^2)-e3*(bl^2)-2*e2*bl*b2; 

f3=-a2*(d2^2)-2*d3*d2*al-2*e3*bl*b2- 

e2*  ((b2'^2)+2*b  1  *b3)+d2*(b3*c  1  +b2*c2+b  1  *c3)+d3*(b2*c  1  +b  1  *c2) ; 
f2=-e3*((b2'^2)+2*bl*b3)-2*e2*b2*b3+d2*(b3*c2+b2*c3)+d3*(b3*cl+b2*c2+bl*c3)- 
a3*(d2'^2)-2*d3*d2*a2-al  *(d3^2); 

fl=b3*c3*d2+d3*(b3*c2+b2*c3)-2*e3*b2*b3-e2*(b3''2)-2*d3*d2*a3-a2*(d3''2); 

f0=b3*c3*d3-a3*(d3''2)-e3*(b3^2); 

coef=[f5  f4  f3  f2  fl  fO]; 

ZG=roots(coef); 


tot(flag)=ZG(5,l); 

end 
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C  PROGRAM  HOPF_lNEW.FOR 

C  EVALUATION  OF  HOPF  BIFURCATION  FORMULAS  USING  THE  SUBOFF 
C  SUBMARINE  MODEL 
C 

IMPLICrr  DOUBLE  PRECISION  (A-H,0-Z) 

REALMS  L,IYY,M,YPDOT,YVDOT,NDl,YRDOT 
REAL*8  YP,YV,YR,NPDOT,NVDOT,NRDOT,NP,NV,NR 
REAL*8  GAMA,U,KVDOT,KRDOT,KPDOT,D 
REAL*8E0,E1,E2,E3,E4,XG,ZG,KR,KP,KV,XG1,ZG1 
C 

REAL*8  Ml  1,M12,M13,M14,M21,M22,M23 
REAL*8  M24,M3 1  ,M32,M33,M34,M41  ,M42,M43,M44 
REAL*8  N1 1,N12,N13,N14,N21,N22,N23,N24 
REAL*8  N3 1  ,N32,N33,N34,N41  ,N42,N43,N44 
REAL*8  LI  1,L12,L13,L14,L21,L22,L23,L24,L31 
REAL*8L32,L33,L34,L15,L16,L17,L25,L26,L27,L35 
REAL*8  L36,L37,L1A,L2A,L3A,L4A,L5A,L6A,L7A,L8A,L9A 
REAL*8  L10A,L1 1  A,L12A,L1,L2,L3,L4,L5,L6,L7 
REAL*8  L8,L9,R1 1,R12,R13,R14,R21,R22,R23,R24 
REAL*8  P11,P12,P13,P21,P22,P23 
DOUBLE  PRECISION  THETA 
C 

DIMENSION  F(4,4),T(4, 4), TINV(4, 4), FV1(4),IV1(4),YYY(4,4) 

DIMENSION  WR(4),WI(4),TSAVE(4,4),TLUD(4, 4), IVLUD(4),SVLUD(4) 
DIMENSION  ASAVE(4,4),A2(4,4),XL(18),HT(18),ZGG(202),FF(4, 4) 
DIMENSION  VEC0(18),VEC1(18),VEC2(18),VEC3(18),VEC4(18),XGG(202) 

C 

INTEGER  I,J,K 
C 

OPEN  (20,FILE=’HOPF25.RES’) 

OPEN  (2 1  ,FILE=DAT25.DAT’,ST  ATUS='OLD') 

C  OPEN  (23,FILE=’HOPFl.RES',STATUS=’OLD') 

C 

WEIGHT=12000.0 
IXX  =1760.0 
lYY  =9450.0 
IZZ  =10700.0 
JXZ  =0.0 
■  rxY  =0.0 
lYZ  =0.0 
L  =17.425 
RHO  =1.94 
CD  =0.5 

C  CD  =0.5*CD*RHO  ??? 

G  =32.2  ' 
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XB  =0.0 
ZB  =0.0 
YG  =0.0 
YB  =0.0 
YDELTAR=0.0 
DELTAR=0.0 
NDELTAR=0.0 
NPROP=0.0 
M  =WEIGHT/G 
B  =  WEIGHT 
W=WEIGHT 

WRITE  (V)’ ENTER  U' 

READ  (*,*.)  U 

WRITE  (*,*)  ’  ENTER  THETA  (DEGREES)' 
READ  (* ,*)  THETA 
C  U  =5.0 

NDl  =  0.5*RHO*L**2 
C  THETA=5 

THETA=THETA*PI/1 80 
OMEGA=U*DTAN(THETA) 

C 

C  DEFINE  HYDRODYNMIIC  COEFnCIENTS 
YPDOT=l  .270E-04*NDl'*L**2 
YVDOT=-5.550E-02*ND1  *L 
YRDOT=l  .240E-03*ND1  *L**2 
YP=3.055E-03*ND1*L 
YPOMEGA=O.DO 

YP=YP+YPOMEGA*OMEGA+M*OMEGA 

YV=-9.310E-02*ND1 

YVOMEGA=O.DO 

yv=yv+yvomega*omega 

YR=-5.940E-02*ND1  *L, 

YROMEGA=O.DO 

YR=YR+YROMEGA*OMEGA 

C 

NPDOT=-3.370E-05*ND  1  *L**3 
NVDOT= 1 .240E-03  *ND  1  *L*  *2 
NRDOT=-3.400E-03*ND1  *L**3 
NP=-8.405E-04*ND  1  *L**2 
NPOMEGA=O.DO 
NP=NP+NPOMEGA*OMEGA 
NV=-1.484E-02*ND1*L 
NVOMEGA=O.DO 

nv=nv+nvomega*omega 

NR=- 1 .640E-02*ND  1  *L**2 
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NROMEGA=O.DO 

NR=NR+NROMEGA*OMEGA 

C 

KPDOT=-1.01E-03*ND1*L**3 
KVD0T=1 .27E-04*ND1*L**2 
KRD0T=-3 .37E-05*ND  1  *L*  *3 
KP=-1.10E-02*ND1*L**2 
KPOMEGA=O.DO 
KP=KP+KPOMEGA*OMEGA 
KV=3.055E-03*ND1*L 
KVOMEGA=O.DO 
KV=KV+KVOMEGA*OMEGA 
.  KR=-8.41E-04*ND1*L**2 
KROMEGA=O.DO 
KR=KR+KROMEGA*OMEGA 
C 

C  DEFINE  THE  DENGTH  X  AND  HEIGHT  H  TERMS  FOR  THE  INTEGRATION 
C  ALLINFEET. 

XL(  1)=-105.9/12.0 
XL(  2)=-104.3/12.0 
XL(  3)=-99.3/12.0 
XL(  4)=-94.3/12.0 
XL(  5)=-87.3/12.0 
XU  6)=-76.8/12.0 
XL(  7)=-66.3/12.0 
XL(  8)=-55.8/12.0 
XL(  9)=72.7/12.0 
XL(10)=79.2/12.0 
XL(11)=83.2/12.0 
XL(12)=87.2/12.0 
XL(13)=91.2/12.0 
XL(14)=95.2/12.0 
XL(15)=99.2/12.0 
XL(16)=101.2/12.0 
XL(I7)=102. 1/12.0 
XL(18)=103.2/12.0 
C 

HT(  1)=  0.000 
HT(  2)=  2.28/12:0 
HT(  3)=  8.24/12.0 
HT(  4)=  13.96/12.0 
HT(  5)=  19.76/12.0 
HT(6)=  25.1/12.0 
HT(  7)=  29.36/12.0 
HT(8)=  31.85/12.0 
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HT(  9)=  31.85/12.0 
HT(10)=  30.00/12.0 
HT(11)=  27.84/12.0 
HT(12)=  25.12/12.0 
HT(13)=  21.44/12.0 
HT(14)=  17.12/12.0 
HT(15)=  12.0/12.0 
HT(16)=  9.12/12 
HT(17)=  6.72/12 
HT(18)=  0.00 
C 

DO  104  K=  1,18 
VEC0(K)=HT(K) 

VEC1(K)=XL(K)*HT(K) 
VEC2(K)=XL(K)*X^K)*HT(K) 
VEC3(K)=XL(K)*XL(K)*XL(K)*HT(K) 
VEC4(K)=XL(K)*XL(K)*XL(K)*XL(K)*HT(K) 
104  CONTINUE 
CALL  TRAP(18,VEC0,XL,E0) 

CALL  TRAP(18,VEC1,XL,E1) 

CALL  TRAP(  1 8,VEC2,XL,E2) 

CALL  TRAP(  1 8,VEC3,XL,E3) 

CALL  TRAP(  1 8,VEC4,XL,E4) 

C 

GAMA=0.001 


C  READ  THE  CRITICAL  VALUES  FOR  XG  AND  ZG  FROM  FILE  DATAl  .DAT 
C 

XGG(1)=0.0 
ZGG(1)=0.016358083 
DO  1 1=2,202 
READ  (21,*)XG,ZG 
C  WRITE(*,*)XG,ZG 
XGG(I)=XG 
ZGG(D=ZG 


C  DETERMINE  [F]  COEFHCIENTS 
C 

D=((M-YVDOT)*(IZZ-NRDOT)*(IXX-KPDOT)) 

&  -((M-YVD0T)*(IXZ-KRD0T)*(-1XZ-NPD0T)) 
&-((M*XG-NVDOT)*(M*XG-YRDOT)*(IXX-KPDOT)) 
&+((M*XG-NVDOT)*(DCZ-KRDOT)*(-M*ZG-YPDOT)) 
&+((-M*ZG-KVD0T)*(M*XG-YRD0T)*(-IXZ-NPD6T)) 
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&  -((-M*ZG-KVDOT)*(IZZ-NRDOT)*(-M*ZG-YPDOT)) 

C 

A11=((IZZ-NRD0T)*(IXX-KPD0T))-((IXZ-KRD0T)*(-IXZ-NPD0T)) 

A12=((-M*XG+YRD0T)*(IXX-KPD0T))+((IXZ-KRD0T)*(-M*ZG-YPD0T)) 

A 1 3=((M*XG-YRD0T)*(-IXZ-NPD0T))-((IZZ-NRD0T)*(-M*ZG-YPD0T)) 
A21=((-M*XG+NVDOT)*(IXX-KPDOT))+((-M*ZG-KVDOT)*(-IXZ-NPDOT)) 
A22=((M-YVDOT)*(IXX-KPDOT))-((-M*ZG-KVDOT)*(-M*ZG-YPDOT)) 
A23=((-M+YVDOT)*(-IXZ-NPDOT))+((M*XG-NVDOT)*(-M*ZG-YPDOT)) 
A31=((M*XG-NVDOT)*(IXZ-KRDOT))-((-M*ZG-KVDOT)*(IZZ-NRDOT)) 
A32=((-M+YVDOT)*(IXZ-KRDOT))+((-M*ZG-KVDOT)*(M*XG-YRDOT)) 
A33=((M-YVDOT)*(IZZ-NRDOT))-((M*XG-NVDOT)*(M*XG-YRDOT)) 
C=================================================== 

C  EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS  • 

C 

F(1,1)=(A11*YV*U+A12*NV*U+A13*KV*U)/D 

F(1,2)=(AI1*(YR*U-M*U)+A12*(-M*XG*U+NR*U)+A13*(M*ZG*U+KR*U))/D 

F(1,3)=(A11*YP*U+A12*NP*U+A13*(KP*U-M*ZG*0MEGA))/D 

F(l,4)=(All*(W-B)*DCOS(THETA)+A12*XG*(W-XB)*B*DCOS(THETA)+ 

&A13*(-ZG*W+ZB)*B*DCOS(THETA))/D 

F(2,1)=(A21*YV*U+A22*NV*U+A23*KV*U)/D 

F(2,2)=(A21*(YR*U-M*U)+A22*(-M*XG*U+NR*U)+A23*(M*ZG*U+KR*U))/D 
F(2,3)=(A2 1  *  YP*U+A22*NP*U+A23*(KP*U-M*ZG*0MEGA))/D 
F(2,4)=(A2 1  *(W-B)*DCOS(THETA)+A22*(XG*W-XB*B)*DCOS(THETA)+ 
&A23*(-ZG*W+ZB*B)*DCOS(THETA))/D 
F(3, 1)=(A3 1  *YV*U+A32*NV*U+A33*KV*U)/D 

F(3,2)=(A31*(YR*U-M*U)+A32*(-M*XG*U+NR*U)+A33*(M*ZG*U+KR*U))/D 
F(3,3)=(A3 1  *  YP*U+A32*NP*U+A33*(KP*U-M*ZG*OMEGA))/D 
F(3,4)=(A3 1  *(W-B)*DCOS(THETA)+A32*(XG*W-XB*B)*DCOS(THETA)+ 
&A33*(-ZG*W+ZB*B)*DCOS(THETA))/D 
F(4,l)=0.0 
F(4,2)=0.0 
F(4,3)=1.0 
F(4,4)=0.0 
C 

D0  11K=1,4 
DO  121=1,4 
ASAVE(K,J)=F(K,J) 

12  CONTINUE 
11  CONTINUE 

CALL  RG(4,4,F,WR,WL1,YYY,IV1,FV1,IERR) 

WR]TE(23, 1007)WR(  1),WR(2),WR(3),WR(4) 
CALLDSOMEG(IEV,WR,WI,OMEGA,CHECK) 

C  WRITE  (*,*)  CHECK 

C  WRITE  (*,*)  WR(2) 
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u  u  u  u  u  u 


C  WRITE  (*  ,*)  WI(4) 

OMEGAO=OMEGA 
DOS  J=l,4 
T(J,1)=  YYY(J,IEV) 

T(J,2)=-YYY(J,IEV+1) 

5  CONTINUE 
IF(IEV.EQ.1.0)GOTO  13 
IF  (IEV.EQ.2.0)  GO  TO  18 
IF  (IEV.EQ.3.0)  GO  TO  14 

C  STOP  3004 

14  D0  6J=1,4 
T(J,3)=YYY(J,1) 

T(J,4)=YYY(J,2) 

6  CONTINUE 
GO  TO  17 

18  D019J=1,4 
T(J,3)=YYY(J,1) 

T(J,4)=YYY(J,4) 

19  CONTINUE 
GO  TO  17 

13  D0  16J=1,4 

T(J,3)=YYY(J,3) 

T(J,4)=YYY(J,4) 

16  CONTINUE 

17  CONTINUE 

NORMALIZATION  OF  THE  CRmCAL  EIGENVECTOR 
CALL  NORM AL(T) 

INVERT  TRANSFORMATION  MATRIX 

D0  2K=1,4 
D0  3  J=l,4 
TINV(K,J)=0.0 
TSAVE(K,J)=T(K,J) 

3  CONTINUE 
2  CONTINUE 

CALLDLUD(4,4,TSAVE,4,TLUD,IVLUD) 

C  D0  4J=1,4 

C  IF  (IVLUD(J).EQ.O)  STOP  3003 

C  4  CONTINUE 

CALL  DILU(4,4,TLUD,IVLUD,SVLUD) 

DO  8  K=l,4 
D0  9J=1,4 
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TINV(K,J)=TLUD(K,J) 

9  CONTINUE 
8  CONTINUE 
C 

C  CHECK  Inv(T)*A*T 
C 

CALL  MULT(TINV,  AS  AVE,T,  A2) 
C 

P1=A2(3,3) 

P2=A2(4,4) 

PEIG1=P1 

PEIG2=P2 

C 

C  DEFINITION  OF  Nij 
C 

N11=TINV(1,1) 

N12=TINV(1,2) 

N13=TINV(1,3) 

N14=TINV(1,4) 

N21=TINV(2,1) 

N22=TINV(2,2) 

N23=TINV(2,3) 

N24=TINV(2,4) 

N31=TINV(3,1) 

N32=TINV(3,2) 

N33=TINV(3,3) 

N34=TINV(3,4) 

N41=TINV(4,1) 

N42=TINV(4,2) 

N43=TINV(4,3) 

N44=TINV(4,4) 

C 

C  DEFINITION  OF  Mij 

C 

M11=T(1.1) 

M12=T(1,2) 

M13=T(1,3) 

M14=T(1,4) 

.  M21=T(2,1) 

M22=T(2,2) 

M23=T(2,3) 

M24=T(2,4) 

M31=T(3,1) 

M32=T(3,2) 

M33=T(3,3) 
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M34=T(3,4) 

M41=T(4,1) 

M42=T(4,2) 

M43=T(4,3) 

M44=T(4,4) 

C 

C  DEFINITION  OF  Lij 
C 

L1=YG*((M3 1  *  *2)+(M2 1  **2)) 

L2=IXY*(M31  **2)+IYZ*M31  *M21+YG*M21  *M1 1 
L3=IXY*M3 1  *M2  l+rYZ*(M2 1  **2)+M*  YG*M  1 1  *M3 1+YG*W*(M4 1  **2)- 
YB*B*(M41** 

&2) 

L4=YG*(2.0*M3 1  *M32+2.0*M21  *M22) 

L5=2.0*IXY*M3 1  *M32+rYZ*(M3 1  *M22+M32*M2 1  )+YG*(M2 1  *M  1 2+M22*M  1 1 ) 

L6=IXY*(M3  1  *M22+M32*M2 1  )+2.0*IYZ*M2 1  *M22+M*  YG*(M  1 1  *M32+M  1 2*M3 1 

)+2. 

&0*(YG*W-YB*B)*M41  *M42 
L7=YG*((M32**2)+(M22**2)) 
L8=IXY*(M32**2)+IYZ*M32*M22+YG*M22*M12 
L9=IXY*M32*M22+IYZ*(M22**2)+M*YG*M12*M32+(YG*W- 
YB*B)*(M42**2) 

C 

C 

L15=(A11/D)*L1+(A12/D)*L2-(A13/D)*L3 
L16=(A1 1/D)*L4+(A12/D)*L5-(A13/D)*L6 
L17=(A11/D)*L7+(A12/D)*L8-(A13/D)*L9 
L25=(A21/D)*L1+(A22/D)*L2-(A23/D)*L3 
L26=(A21/D)*L4+(A22/D)*L5-(A23/D)*L6 
L27=(A21/D)*L7+(A22/D)*L8-(A23/D)*L9 
L35=(A3 1/D)*L1+(A32/D)*L2-(A33/D)*L3 
L36=(A3 1/D)*L4+(A32/D)*L5-(A33/D)*L6 
L37=(A3 1/D)*L7+(A32/D)*L8-(A33/D)*L9 
C 

C=CD/(6.0*GAMA) 

L1A=C*(E0*(M1 1**3)+3.0*E1*(M1 1**2)*M21+3.0*E2*M1 1*(M21**2)+E3*(M21 
&**3))-((l-0/6.0)*(W-B)*DCOS(THETA))*(M41**3) 

L2A=C*(E1*(M11**3)+3.0*E2*(M11**2)*M21+3.0*E3*M11*(M21**2)+E4*(M21 
&**3))-((l  .0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M41  **3) 

L3  A=((  1 .0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M4 1  **3) 
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L4A=C*(3.0*E0*(Mll**2)*M12+3.0*El*((Mll**2)*M22+2*Mll*M12*M21)+3.0 
&*E2*(M12*(M21**2)+2.0*M21*M22*Mll)+3.0*E3*(M21**2)*M22)-((1.0/6.0) 
&*(W-B)*DCOS(THETA))*3.0*(M41  **2)*M42 

L5A=C*(3.0*El*(Mll**2)*M12+3.0*E2*((Mll**2)*M22+2*Mll*M12*M21)+3.0 

&*E3*(M12*(M21**2)+2*M21*M22*M11)+3.0*E4*(M21**2)*M22)-((1.0/6.0)*(X 

&G*W-XB*B)*DCOS(THETA))*3.0*(M41**2)*M42 

L6A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*(M41**2)*M42 

L7A=C*(3*E0*M11*(M12**2)+3.0*E1*((M12**2)*M21+2*M11*M12*M22)+3.0*E 
&2*((M22**2)*M11+2.0*M21*M22*M12)+3.0*E3*M21*(M22**2)H(1.0/6.0)*(W 
&-B)*DCOS(THETA))*3.0*M4 1  *(M42**2) 

L8A=C*(3.0*El*Mll*(M12**2)+3.0*E2*((M12**2)*M21+2.0*Mll*M12*M22)+3 

&.0*E3*((M22**2)*M11+2*M21*M22*M12)+3.0*E4*M21*(M22**2)) 

L9A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*3.0*M41*(M42**2) 

L10A=-C*(E0*(M12**3)+3.0*E1*(M11**2)*M21+3.0*E2*M12*(M22**2)+E3*(M 

&22**3))-((1.0/6.0)*(W-B)*DCOS(THETA))*(M42**3) 

LI  1  A=-C*(E1*(M12**3)+3.0*E2*(M1 1**2)*M21+3.0*E3*M12*(M22**2)+E4*(M 

&22**3))-((1.0/6.0)*(XG*W-XB*B)*DCOS(THETA))*(M42**3) 

L12A=((1.0/6.0)*(ZG*W-ZB*B)*DCOS(THETA))*(M42**3) 

C 

LI  1=(-A1 1/D)*L1  A+(A12/D)*L2A+(A13/D)*L3A 
L12=(-A11/D)*L4A+(A12/D)*L5A+(A13/D)*L6A 
L13=(-A11/D)*L7A+(A12/D)*L8A+(A13/D)*L9A 
L14=(-A1 1/D)*L10A+(A12/D)*L1 1  A+(A13/D)*L12A 
C 

L21=(-A21/D)*L1A4-(A22/D)*L2A+(A23/D)*L3A 

L22=(-A21/D)*L4A+(A22/D)*L5A+(A23/D)*L6A 

L23=(-A21/D)*L7A+(A22/D)*L8A+(A23/D)*L9A 

L24=(-A21/D)*L10A+(A22/D)*L11A+(A23/D)*L12A 

C 

L31=(-A31/D)*L1A+(A32/D)*L2A+(A33/D)*L3A 

L32=(-A31/D)*L4A+(A32/D)*L5A+(A33/D)*L6A 

L33=(-A31/D)*L7A+(A32/D)*L8A+(A33/D)*L9A 

L34=(-A31/D)*L10A+(A32/D)*L11A+(A33/D)*L12A 

C 

R1 1=(N1 1*L1 1)+(N12*L21)+(N13*L31) 

R12=(N11*L12)+(N12*L22)+(N13*L32) 

R13=(N1 1*L13)+(N12*L23)+(N13*L33) 

R14=(N1 1*L14)+(N12*L24)+(N13*L34) 

R21=(N21*L1 1)+(N22*L21)+(N23*L3 1) 

R22=(N2 1  *L1 2)+(N22*L22)+(N23  *L32) 

R23=(N2 1  *L1 3)+(N22*L23)+(N23*L33) 
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non  o'  on  o  oo  non 


R24=(N2 1  *L14)+(N22*L24)+(N23*L34) 

C 

P11=(N11*L15)+(N12*L25)+(N13*L35) 
P12=(N11*L16)+(N12*L26)+(N13*L36) 
P13=(N11*L17)+(N12*L27)+(N13*L37) 
P21=(N21*L15)+(N22*L25)+(N23*L35) 

P22=(N2 1  *L16)+(N22*L26)+(N23*L36) 
P23=(N21*L17)+(N22*L27)+(N23*L37) 

EVALUATE  DALPHA  AND  DOMEGA 

ZGR=ZGG(I) 

ZGL=ZGG(I-1) 

ZGl  =ZGR 
XGl  =XGG(I) 

CALL  FMATRIX(ZG1  ,XG1  ,FF,U, THETA) 

CALL  RG(4,4,FF,WR,WI,0,  YYY,IV  1  ,FV  1  ,ffiRR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

omEgr=freq 

ZGl  =ZGL 
XGl  =XGG(I-1) 

CALL  FMATRIX(ZG1  ,XG1  ,FF,U, THETA) 

CALL  RG(4,4,FF,WR,WI,0,YYY,IV1  ,FV  1  ,IERR) 

CALL  DSTABL(DE0S,WR,WI,FREQ) 

ALPHL=DEOS 
OMEGL=FREQ 

DALPHA=(ALPHR-ALPHL)/(ZGR-ZGL) 
DOMEGA=(OMEGR-OMEGL)/(ZGR-ZGL) 

EVALUATION  OF  HOPF  BIFURCATION  COEFnCIENTS 

.  COEF1=(1.0/8.0)*(3.0*R11+R13+R22+3.0*R24) 
COEF2=(1.0/8.0)*(3.0*R11+R23-R12-3.0*R14) 

C  AMU2  =-COfiFl/(8.0*DALPHA)  ???????? 

C  BETA2=0.25*COEF1  ??????? 

C  TAU2  =-(COEF2-DOMEGA*COEFl/DALPHA)/(8.0*OMEGAO) 
C  PER  =2.0*3.1415927/OMEGA0 
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C  PER  =PER*U/L 

WRITE  (20,2001)XG,ZG,COEF1  ,DALPHA,OMEGAO,PEIG1  ,PEIG2 
C  WRITE  (20,2001)COEF1 
1  CONTINUE 
C 

STOP 

2001  FORMAT  (7E14.5) 

4001  FORMAT  (3F15.5) 

1007  FORMAT  (4E14.5) 

END 
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